library("openxlsx")
G1=read.xlsx("D:/浏览器下载/data for class3/newfile3.xlsx")
G1[G1==0]=NA
row.names(G1)=G1$Row.names
G1=G1[,-1]
G2=log2(G1)
colnames(G1)
group=gsub("\\d","",colnames(G1))#提取分组信息
table(group)
P=names(as.data.frame(G2))#创建新的表格用来存接下来的数据
#利用ANOVA检验分别求所有蛋白在3组之间的P值
for (i in 1:100) {
  pro=unlist(G2[i,])#取第i行的蛋白质
  #提取各组数据
  posctl=group=="ctl"
  ctl=pro[posctl]
  posad=group=="ad"
  ad=pro[posad]
  posasym=group=="asym"
  asym=pro[posasym]
  #统计每组中非0的个数
  posx=!is.na(ctl)
  cltnum=sum(posx)
  posx2=!is.na(ad)
  adnum=sum(posx2)
  posx3=!is.na(asym)
  asymnum=sum(posx3)
  {if(cltnum> 3& adnum>3 & asymnum>3){
    result=oneway.test(pro~group)#ANOVA检验
    #把输出的结果作为新的一列保存
    P[i]=result$p.value
  }
  else{P[i]=NA}}
  
  
}
print(P)
#GO富集分析
library("BiocManager")
library("scatterpie")
BiocManager::install("clusterProfiler")

library("clusterProfiler")
library("org.Hs.eg.db")
library("ggplot2")
load("D:/浏览器下载/class5_volcano.RData")
#筛选p<0.05的蛋白
DEG=prostat$ID[prostat$P<0.05] 
allID=prostat$ID[!is.na(prostat$ID)]
#基因背景为列表基因，对其功能BP富集分析
enall= enrichGO(DEG,OrgDb = "org.Hs.eg.db",keyType = "SYMBOL",
                   ont = "BP",pvalueCutoff = 1, qvalueCutoff = 1, universe=allID)
#条形图展示
barplot(enall,showCategory = 20)
#气泡图展示
dotplot(enall,showCategory = 20)
